

xtset id year               

/************************************Table 1. Basic Statistics.*****************************************************************************/


tabstat TCE  CEP   IPRP ES TP RD PS  Trade IS  Urb   EG  , stats(count mean sd min max) c(s) f(%6.2f)
logout, save(mytable) word replace: tabstat TCE  CEP   IPRP ES TP RD PS  Trade IS  Urb   EG, stats(count mean sd min max) c(s) f(%6.2f)

/************************************Figure 2. Kernel density estimation plot.*****************************************************************************/

kdensity TCE if year==2008 , addplot( (kdensity  TCE if year==2009,clp(longdash)) (kdensity  TCE if year==2010,clp(dash_dot))  (kdensity  TCE if year==2011, clp(shortdash))(kdensity  TCE if year==2012,clp(dash_dot)) (kdensity  TCE if year==2013,clp(dash_dot)) (kdensity  TCE if year==2014,clp(dash_dot)) (kdensity  TCE if year==2015,clp(dash_dot)) (kdensity  TCE if year==2016,clp(dash_dot)) (kdensity  TCE if year==2017,clp(dash_dot)) (kdensity  TCE if year==2018,clp(dash_dot)) (kdensity  TCE if year==2019,clp(dash_dot)) (kdensity  TCE if year==2020,clp(dash_dot)) ) legend( label (1 "2008 ") label (2 "2009 ") label (3 "2010 ") label (4 "2011 ") label (5 "2012 ") label (6 "2013 ") label (7 "2014 ") label (8 "2015 ") label (9 "2016 ") label (10 "2017 ") label (11 "2018 ") label (12 "2019 ") label (13 "2020 ") size(small) position(1) row(4) ring(0) subtitle(year)  region(fcolor(gs15)))  title("Kernel density estimate")  xtitle("TCE") 
 
kdensity CEP if year==2008 , addplot( (kdensity  CEP if year==2009,clp(longdash)) (kdensity  CEP if year==2010,clp(dash_dot))  (kdensity  CEP if year==2011, clp(shortdash))(kdensity  CEP if year==2012,clp(dash_dot)) (kdensity  CEP if year==2013,clp(dash_dot)) (kdensity  CEP if year==2014,clp(dash_dot)) (kdensity  CEP if year==2015,clp(dash_dot)) (kdensity  CEP if year==2016,clp(dash_dot)) (kdensity  CEP if year==2017,clp(dash_dot)) (kdensity  CEP if year==2018,clp(dash_dot)) (kdensity  CEP if year==2019,clp(dash_dot)) (kdensity  CEP if year==2020,clp(dash_dot)) ) legend( label (1 "2008 ") label (2 "2009 ") label (3 "2010 ") label (4 "2011 ") label (5 "2012 ") label (6 "2013 ") label (7 "2014 ") label (8 "2015 ") label (9 "2016 ") label (10 "2017 ") label (11 "2018 ") label (12 "2019 ") label (13 "2020 ") size(small) position(1) row(4) ring(0) subtitle(year)  region(fcolor(gs15)))  title("Kernel density estimate")  xtitle("CEP") 

kdensity IPRP if year==2008 , addplot( (kdensity  IPRP if year==2009,clp(longdash)) (kdensity  IPRP if year==2010,clp(dash_dot))  (kdensity  IPRP if year==2011, clp(shortdash))(kdensity  IPRP if year==2012,clp(dash_dot)) (kdensity  IPRP if year==2013,clp(dash_dot)) (kdensity  IPRP if year==2014,clp(dash_dot)) (kdensity  IPRP if year==2015,clp(dash_dot)) (kdensity  IPRP if year==2016,clp(dash_dot)) (kdensity  IPRP if year==2017,clp(dash_dot)) (kdensity  IPRP if year==2018,clp(dash_dot)) (kdensity  IPRP if year==2019,clp(dash_dot)) (kdensity  IPRP if year==2020,clp(dash_dot)) ) legend( label (1 "2008 ") label (2 "2009 ") label (3 "2010 ") label (4 "2011 ") label (5 "2012 ") label (6 "2013 ") label (7 "2014 ") label (8 "2015 ") label (9 "2016 ") label (10 "2017 ") label (11 "2018 ") label (12 "2019 ") label (13 "2020 ") size(small) position(1) row(4) ring(0) subtitle(year)  region(fcolor(gs15)))  title("Kernel density estimate")  xtitle("IPRP") 


 
/*********************************************Table 2. Baseline results**************************************************************************/
reghdfe  TCE     IPRP            , absorb(id year ) vce(cluster id)   
est store a
reghdfe  CEP    IPRP            , absorb(id year ) vce(cluster id)   
est store b
reghdfe  TCE     IPRP   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store c
reghdfe  CEP    IPRP   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store d
esttab a b c d     using reg1.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap



/***********************VIF******************/

reg TCE     IPRP   Trade IS  Urb   EG  
estat vif
reg CEP    IPRP   Trade IS  Urb   EG 
estat vif

/***********Hausman Test*****************/

xtreg TCE     IPRP   Trade IS  Urb   EG   ,fe 
esti store FE
xtreg TCE     IPRP   Trade IS  Urb   EG   ,re 
esti store RE
hausman FE RE,constant sigmamore

xtreg CEP     IPRP   Trade IS  Urb   EG ,fe
esti store FE
xtreg CEP     IPRP   Trade IS  Urb   EG ,re
esti store RE
hausman FE RE,constant sigmamore

/*********************************************Robustness Test**************************************************************************/
/***********************Table 3. Robustness test results.******************/

reghdfe  TCE     IPRIGCI   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store a
reghdfe  CEP    IPRIGCI   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store b
reghdfe  TCE     IPRP   Trade IS  Urb   EG     if year<=2019  , absorb(id year ) vce(cluster id)   
est store e
reghdfe  CEP    IPRP  Trade IS  Urb   EG    if year<=2019  , absorb(id year ) vce(cluster id)   
est store f
reghdfe  CO2     IPRP   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store c
reghdfe  BTFP     IPRP  Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store d

esttab a b c d e f   using reg2.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap

/***********************Table 4. 2SLS results.******************/
tab year, gen(dyear)

xtivreg2   TCE    Trade IS  Urb   EG  dyear1-dyear13  ( IPRP =  RL  )  ,fe first 
est store a
xtivreg2   CEP   Trade IS  Urb   EG dyear1-dyear13  ( IPRP =  RL  )  ,fe first 
est store b
xtivreg2   TCE    Trade IS  Urb   EG  dyear1-dyear13  ( IPRP =  l.IPRP  )  ,fe first 
est store c
xtivreg2   CEP   Trade IS  Urb   EG dyear1-dyear13  ( IPRP =  l.IPRP  )  ,fe first 
est store d
xtivreg2  TCE    Trade IS  Urb   EG    LnRGDP1 dyear1-dyear13  ( IPRP =  RL l.IPRP  )  ,fe first 
est store e
xtivreg2   CEP   Trade IS  Urb   EG dyear1-dyear13  ( IPRP =  RL l.IPRP  )  ,fe first 
est store f
esttab a b c d e f   using reg3.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap


/*********************************************Table 5. Heterogeneity analysis results.**************************************************************************/
reghdfe   TCE     IPRP   Trade IS  Urb   EG    HIC  c.IPRP#c.HIC       , absorb(id year ) vce(cluster id)   
est store a
reghdfe   TCE     IPRP   Trade IS  Urb   EG    UMIC  c.IPRP#c.UMIC       , absorb(id year ) vce(cluster id)   
est store b
reghdfe   TCE     IPRP   Trade IS  Urb   EG    LMIC  c.IPRP#c.LMIC       , absorb(id year ) vce(cluster id)   
est store c
reghdfe   CEP    IPRP   Trade IS  Urb   EG    HIC  c.IPRP#c.HIC       , absorb(id year ) vce(cluster id)   
est store d
reghdfe   CEP    IPRP   Trade IS  Urb   EG    UMIC  c.IPRP#c.UMIC       , absorb(id year ) vce(cluster id)   
est store e
reghdfe   CEP    IPRP   Trade IS  Urb   EG    LMIC  c.IPRP#c.LMIC       , absorb(id year ) vce(cluster id)   
est store f
esttab a b c d e f   using reg6.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap



/*********************************************Table 6. Mechanism analysis result.**************************************************************************/


reghdfe  ES     IPRP   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store a
reghdfe  TP    IPRP   Trade IS  Urb   EG        , absorb(id year ) vce(cluster id)   
est store b
reghdfe  EG     IPRP   Trade IS  Urb          , absorb(id year ) vce(cluster id)   
est store c
reghdfe  TCE  ES  IPRP   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store d
reghdfe  CEP ES  IPRP   Trade IS  Urb   EG        , absorb(id year ) vce(cluster id)   
est store e
reghdfe  TCE  TP IPRP   Trade IS  Urb   EG      , absorb(id year ) vce(cluster id)   
est store f
reghdfe  CEP TP IPRP   Trade IS  Urb   EG        , absorb(id year ) vce(cluster id)   
est store g
reghdfe  TCE  ES  IPRP   Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store h
reghdfe  CEP ES  IPRP   Trade IS  Urb   EG        , absorb(id year ) vce(cluster id)   
est store i

esttab a b c d e f g h i  using reg4.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap



/*********************************************Table 7. Moderating effect analysis results**************************************************************************/

xtreg   TCE     c_IPRP   Trade IS  Urb   EG c_RQ  c.c_IPRP#c.c_RQ    i.year, fe r     
est store a
xtreg   CEP    c_IPRP   Trade IS  Urb   EG c_RQ  c.c_IPRP#c.c_RQ    i.year, fe r   
est store b

xtreg   TCE     c_IPRP   Trade IS  Urb   EG c_PS  c.c_IPRP#c.c_PS    i.year, fe r     
est store c
xtreg   CEP    c_IPRP   Trade IS  Urb   EG c_PS  c.c_IPRP#c.c_PS    i.year, fe r   
est store d

esttab a b c d   using reg5.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap

/*********************************************4.6.Further analysis**************************************************************************/


 /*********Figure 3. Nonlinear fitting diagram***********/


reg  TCE     IPRP   Trade IS  Urb   EG       
estat ovtest,rhs
reg   CEP    IPRP   Trade IS  Urb   EG      
estat ovtest,rhs


reghdfe  TCE     IPRP IPRP2  Trade IS  Urb   EG      , absorb(id year ) vce(cluster id)   
est store a
reghdfe  CEP     IPRP IPRP2  Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   
est store b
esttab a b    using reg7.rtf, r2 ar2 se star(* 0.1 ** 0.05 *** 0.001) replace nogap

reghdfe  TCE     IPRP IPRP2  Trade IS  Urb   EG       , absorb(id year ) vce(cluster id)   

global tp: dis %3.1f -_b[IPRP]/(2*_b[IPRP2]) //
dis "Turn point = " $tp 

local b0 = _b[_cons]
 . local b1 = _b[IPRP]
 . local b2 = _b[IPRP2]
 . sum IPRP
 . local min:  dis %3.1f r(min)
 . local max:  dis %3.1f r(max)
#delimit
twoway ( function y = `b2'*x^2 + `b1'*x + `b0',range(`min'  `max') ) ( function y = `b2'*x^2 + `b1'*x + `b0', range(`max' 2) lp(dash) ), ytitle("TCE") xtitle("IPRP") xline(`min',lc(red) lw(*1.5))  xline($tp,  lp(dash) lc(green))  xline(`max',lc(red) lw(*1.5))  xlabel(`min' "Min (`min')" 0.7(0.1)2.5 `max' "Max (`max')", angle(30))   legend(off);
  #delimit cr
  graph export "fig_Ushape.png", replace  


/*********************************************Figure 4. Marginal effect diagram**************************************************************************/
reghdfe  TCE   c.IPRP##c.IPRP Trade IS  Urb   EG         , absorb(id year ) vce(cluster id) $xx
sum IPRP
global min: dis %4.1f r(min)
global max: dis %4.1f r(max)
margins, dydx(IPRP) at(IPRP=($min 0.7(0.1)2.5 $max))

marginsplot, xlabel(, format(%3.1f) angle(60)) //
      graph export "$Out\fig_Wage_marginplot_A.png", replace

      marginsplot, addplot(function y=0, range($min $max)  ///
                   lcolor(red) lpattern(dash) legend(off)) ///
                   xlabel(, format(%3.1f) angle(60)) ///
                   ylabel(, format(%2.1f) angle(0))
      graph export "$Out\fig_Wage_marginplot_B.png", replace

/*********************************************end**************************************************************************/